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Abstract. The Generalized Aeroelastic Analysis Method (GAAM) is applied to the analysis 
of three well-studied checkcases: restrained and unrestrained airfoil models, and a wing 
model. An eigenvalue iteration procedure is used for converging upon roots of the complex 
stability matrix. For the airfoil models, exact root loci are given which clearly illustrate the 
nature of the flutter and divergence instabilities. The singularities involved are enumerated, 
including an additional pole at the origin for the unrestrained airfoil case and the emergence 
of an additional pole on the positive real axis at the divergence speed for the restrained airfoil 
case. Inconsistencies and differences among published aeroelastic root loci and the new, 
exact results are discussed and resolved. The generalization of a Doublet Lattice Method 
computer code is described and the code is applied to the calculation of root loci for the wing 
model for incompressible and for subsonic flow conditions. The error introduced in the 
reduction of the singular integral equation underlying the unsteady lifting surface theory to a 
linear algebraic equation is discussed. Acknowledging this inherent error, the solutions of the 
algebraic equation by GAAM are termed ‘exact.’ The singularities of the problem are 
discussed and exponential series approximations used in the evaluation of the kernel function 
shown to introduce a dense collection of poles and zeroes on the negative real axis. Again, 
inconsistencies and differences among published aeroelastic root loci and the new, ‘exact’ 
results are discussed and resolved. In all cases, aeroelastic flutter and divergence speeds and 
frequencies are in good agreement with published results. The GAAM solution procedure 
allows complete control over Mach number, velocity, density, and complex frequency. Thus 
all points on the computed root loci can be matched-point, consistent solutions without 
recourse to complex mode tracking logic or dataset interpolation, as in the k and p-k solution 
methods. 

1 INTRODUCTION 

Aeroelastic divergence of a restrained lifting surface or vehicle is a static stability problem 
that can be solved easily and accurately by eigenvalue methods for general, real matrices (e.g. 
pp. 431-440 of Ref. 1). The stability of an unrestrained vehicle, however, is a dynamic 
problem that cannot be solved accurately by static or quasi-steady methods. A number of 
recent publications have brought attention to differences and inconsistencies in calculations of 
the flutter and divergence behaviors of restrained airfoil 2 ' 4 , unrestrained airfoil 3 ’ 5 ’ 6 , and 
cantilevered wing 3 ’ 7 ’ 8 models when analyzed by dynamic stability methods. The use of a 
number of approximations to the unsteady airloads has led to this situation. Also, the various 
procedures of incorporating available, harmonic, unsteady airloads into the aeroelastic 
analyses, the methodologies of ‘root-sorting’ and ‘mode-tracking’ for reporting results, and 
the use of ‘decay rates’ as ‘damping values’ for nonoscillatory, real roots (with attendant 
‘bifurcations’, ‘jumps’, and ‘activation of lag roots’) contribute to the differences. 
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The unrestrained three Degree Of Freedom (DOF) airfoil model shown in Fig. 1 was studied 
by Rodden and Bellinger 5 using the p and p-k methods and is Example HA145A in Ref. 9. 
The case has also been studied by Chen 3 using the g-method and by van Zyl 6 using four forms 
of the p-k method. The restrained 2 DOF airfoil model (obtained by eliminating the 
‘fuselage’ in Fig.l, or by letting the the fuselage mass approach infinity) was studied by 
Rodden and Bellinger 2 using the p, p-k, and k methods. The case has also been studied by 
Chen - ’ and by van Zyl 4 . Both cases have been studied for two locations of the center of 
gravity. For the forward location, the divergence speed, U D , is lower than the flutter speed, 
U f . Also, for the 3 DOF unrestrained case, a low frequency oscillatory mode instability 

appears instead of the well-known quasi-static divergence behavior. Since the mode “finds its 
origin in a tendency to static divergence,” 5 it was termed ‘dynamic divergence’. Small 
differences in the flutter and divergence speeds are reported for these cases. More 
significantly, differences are reported in the composition, origin, and continuity of the system 
modes for speeds well below and above U D and U f , and for highly damped or undamped 

modes. Thus, there are discussions of speeds at which aerodynamic lag roots ‘become 
active’ 3,10 , speeds at which two damped, oscillatory modes merge and are no longer found for 
higher velocities 4 , and differences over the origin of the divergence root 2,4 ' 6 (i.e., whether it 
derives continuously from a structural mode or from an aerodynamic lag root). 



^ rf-ii , « „ , ...... aerodynamic strip modeling. 

Figure 1. Three degree of freedom airfoil and J 

‘fuselage.’ 


The jet transport wing model, the BAH wing, shown in Fig. 2 was introduced by 
Bisplinghoff, Ashley, and Halfman 1 and has become a standard checkcase for flutter and 
divergence analysis of wings. It is Example HA145B in Ref. 9. Rodden and Stahl 7 studied 
the BAH wing using the transient (p) method with airloads from strip theory and in Ref. 9 it is 
studied using the p-k method with the loads computed by the Doublet Lattice Method 11 
(DLM). The wing has also been studied by Chen 3 using the g-method and by van Zyl s using 
three forms of the p-k method. There is agreement for the general structure of the 
instabilities: a bending-torsion flutter mode and a real divergence root. Again, small 
differences in U D and U f are noted, depending upon details of the calculations. Also, the 

discussions and differences noted above for the airfoil cases regarding the composition, 
origin, and continuity of the system modes are also present for the BAH wing case. 

All of the above calculations have been based upon unsteady aerodynamic theories which 
assume purely harmonic structural oscillations. This leads directly to many of the issues 
underlying the results discussed above. A generalization of the harmonic oscillation 
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assumption has been available and will be shown to resolve the issues and differences. 
Edwards 12 14 and Edwards, Ashley, and Breakwell 15 proved the validity of the unsteady, 
linear, potential equation for arbitrary complex values of the Laplace transform variable, s, 
giving examples for Two-Dimensional (2-D) incompressible 12 ’ 13 ’ 15 , subsonic 14 , and 
supersonic 12,15 speeds and for Three-Dimensional (3-D) subsonic 14 speeds. Examples of exact 
aeroelastic root loci for general complex values of s illustrating flutter 12 ’ 13 ’ 15 and 
divergence 12 ’ 15 instabilities are given for 2-D incompressible and supersonic flow. For 
incompressible flow, it was shown that, above the divergence speed, an additional pole (the 
divergence mode) is found on the positive real axis. The mode exists in addition to those 
deriving from the structural degrees of freedom. Edwards 14 also showed that existing 
computer codes, developed for the calculation of harmonic oscillation airloads, could be 
generalized in a straightforward manner. A DLM computer code was generalized, and 
damped and undamped airloads were calculated for a swept and tapered wing. Cunningham 
and Desmarais 16 and Cunningham 17 performed similar generalizations of 2-D and 3-D 
unsteady potential theory Kernel Function Method codes for subsonic and supersonic flow, 
respectively. They present generalized airloads for an oscillating airfoil and generalized 
flutter analyses for several wing configurations. 

The purpose of this paper is to present new analyses of the restrained and unrestrained airfoil 
models, and the BAH wing model using a Generalized Aeroelastic Analysis Method. 
Generalized unsteady airloads (incompressible and subsonic) are used to calculate the exact 
root loci of these models throughout the complex s-plane, including their flutter and 
divergence behaviors. Along with attention to the number and nature of the singularities 
involved, the new analysis gives a complete, intuitive explanation of the stability of the 
models and new insights into the capabilities and limitations of the computational methods in 
use. 

2 AEROELASTIC SOLUTION METHODS 

Aeroelastic stability analyses are based upon the Laplace transformed flutter equation 

[M S 5 2 + J B^ + K S -^(T,M)]{^(5)} = 0 (1) 

where A/., B s , and K s are the structural Mass, Damping, and Stiffness matrices and A is the 
Generalized Aerodynamic Force (GAF) matrix which is a function of the generalized reduced 
frequency, J = sb/U and Mach number, M . Hassig ls summarizes aeroelastic solution 
methods which have become well known: the p method (‘transient’ method 2 ’ 7 ), the k method 
(‘American’ method, ‘U-g’ method, ‘K’ and ‘KE’ methods in MSC/NASTRAN 9 ), and the p-k 
method (‘British’ method, ‘PK’ method in MSC/NASTRAN 9 ). All of these methods have 
been based upon the assumption that the required unsteady aerodynamic loads, A, are 
available, except in some special cases, only for purely harmonic oscillation of the structural 
degrees of freedom. References 12-15 show that this has been an unnecessary restriction and 
that airloads for arbitrary complex values of frequency, p (or.y = a+/<y herein), can be 
computed by suitably generalized computational methods and codes. 

Returning to the development of solution methods, the p-k method has been further 
differentiated by the treatment of aerodynamic damping. Rodden, Harder, and Bellinger 19 
introduced the A matrix as 
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A=(\pcUQ I lk)p + \pU 2 Q R 


( 2 ) 


where Q = Q R +iQ‘ is the matrix of GAFs for purely harmonic oscillations at a given reduced 
frequency ~s =ik = iO)b/U . In his ‘g-method’, Chen 3 introduces the A matrix as 


A = ±pU 2 Q'(ik)g+±pcU 2 Q(ik) 
Finally, van Zyl 8 uses a variation on the g-method, 


A=-'pUb 


mm 

dk 


p+lpU 2 


Q(ik)-k d( @m 


(3) 

(4) 


with special attention to numerical problems at k = 0 by means of spline approximation of the 
elements of Q and its derivative. The use of derivatives of the GAF matrix in Eqs. 3 and 4 
makes these methods questionable for use in divergence analysis, with k = 0. Subsonic 
airloads have a logarithmic branch point at the origin where their derivatives are undefined. 


All of the checkcases to be considered are for incompressible flow, M = 0 . For the BAH 
wing model, the generalized DLM code described below is also directly applicable for 
subsonic, compressible flow and example calculations are given. For the 2-D airfoil cases the 
Q matrix loads are available in closed form (e.g., Eqs. 5-3 1 1,3 12 of Ref. 1). They involve the 
Theodorsen function, C(ik), which can be utilized in tabulated form (for a series of reduced 
frequencies, k, and with interpolation for intermediate values) or in approximate forms, such 
as that of W. P. Jones 20 wherein the Wagner function is approximated by a two term 
exponential series curve fit 


*(Ut/b) » 1- exp (-P n Wb) (5) 

C(ik) and Q>(Ut/b) were known to be a Fourier transform pair, enabling the approximation 
of Eq. 5 to be incorporated into the flutter equation (1) as two additional linear, constant 
coefficient. Tag’ equations. These properties of the lag equations suggested that the 
approximation could be used for arbitrary complex values of reduced frequency, which is the 
essence of the p method for flutter analysis. 

As summarized by Rodden and Johnson 9 these methods and approximations lead to forms of 
Eq. 1 which can be solved by linear matrix eigenvalue methods for general, real (Eq. 2) and 
complex (Eqs. 3 and 4) matrices. In the k-method, the introduction of structural damping, g, 
leads to an NxN matrix equation for the ‘eigenvalue’ variable 

p 2 =(-U 2 /(l + ig)) = a n +ib n ,n = l,...,N. Velocity and damping values are recovered as 

U n = + b n 2 ) / a n and g n = -b n / a n . Interpretation of the results leads to complex mode 

tracking and root sorting logic and solution iterations are required in order to obtain consistent 
‘matched point’ flutter solutions. The p-k method (p method) leads to a 2Nx2N (2Nx2N + 2) 
matrix equation for the variable p = co(g + i) where the viscous damping is g =jg. Real 

roots resulting from the lag equations in the p method, and the occurrence of real roots for the 
p-k method, have led to the use 9 of the decay rate of real, non-oscillatory roots as damping 
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values in which g is assigned the value 2Re(/>)Z>/£/(ln2) . This facilitates some continuity 

in plots of the frequencies and dampings; however, many details of the resulting 
U — g and U-f plots are incorrect as will be shown below. 

Ref. 21 discusses the problems of mode-tracking and mode-switching which surface in these 
aeroelastic solution methods and investigates two methods of improving solution performance 
by using information from the eigenvectors. Ref. 22 introduces a ‘piecewise aerodynamic 
flutter method,’ an alternative to the p-k method, in order to address problems with solution 
convergence and the difficult and time-consuming construction of state-space models. 

3 MODEL DESCRIPTIONS 


3.1 Airfoil Models 

The transformed equations of motion for the 3 DOF airfoil of Fig. 1 are presented in Eq.6: 
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The bending and torsion springs are attached to a mass (fuselage), free to move vertically. 
Two examples are considered, differing only in the location of the center of gravity: in 
example 1 it is at 37% chord and in example 2 it is at 45% chord. Both examples assume 
incompressible flow. The aerodynamic center is at 25% chord and the elastic axis is at 40% 
chord. The chord is 6 feet, the nondimensional radius of gyration is >« = 0.25 , and the mass 
ratio // = 20. The static un-balances are x a = -0.06 for example 1 and x a = 0. 10 for example 
2. The uncoupled bending and torsion frequencies are 0) h = 10 rad/sec and 0) a =25 rad/sec. 
Equal viscous damping of g =\g = {0.03 is applied to the first two modes. The fuselage 
mass is assumed to be equal to the airfoil mass, m f = m w . The 2 DOF airfoil model is 
obtained from Eq. 6 by eliminating the row and column associated with the fuselage, h f , 
keeping all the airfoil properties unchanged. 

For incompressible flow, the airload coefficients in Eq. 6 are 


r ~c,(s)b 
, cjs) 




J 


\_ M nc s + B nc s + K nc + C(s )R(S 2 s -t-S)) 



(7) 


The matrices in Eq. 7 are defined in Ref. 13. The generalized Theodorsen function is 


C(j) = ¥ ^ ,s 1 • 

K^D + K^s)’ 



( 8 ) 


where Ko and Kj are modified Bessel functions which are defined and analytic throughout the 
s -plane except for logarithmic branch points at the origin and at infinity. A branch cut along 
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the negative real axis makes them single-valued. Combining Eqs. 6 and 7 and defining the 
state vector X T = (h w ,a,h f ,sh w ,sa,sh f ), the equations of motion can be written as 13 


{sI-^(s)}JST(s) = 0 

where 


A(s) = 


0 I 

M~'K(s) M~'B(s ) 


( 9 ) 

( 10 ) 


Note that AT(s)and B(s) are functions of C(s). Thus the stability matrix, A(s) , has non- 
constant, complex elements and attention must be paid to the logarithmic singularity at the 
origin. Thus, common eigenvalue methods of determining system stability using computer 
codes for real, general matrices cannot be applied to Eq. 9. Instead, a shareware family of 
routines for determining eigenvalues and eigenvectors of complex, general matrices was 
utilized, along with an eigenvalue iteration method to converge upon, and to track the 
singularities (e.g., poles) of Eq. 9. Thus, starting from an initial search location, s„, the 

eigenvalues X T = (s l ,...,s 6 } of A(s „) are determined and interrogated for the value of 

s' closest to s n . An iterative relaxation method was used to converge upon a root of Eq. 9: 


Sn + i=Sn+r fac (s‘-s n ) ( 11 ) 

In this manner, exact root loci of the modes of Eq. 9 are calculated by incrementing U 
following the determination of a converged root. Typical starting locations are the coupled, 
wind-off structural frequencies. The method was very robust for the complex structural 
modes; relaxation factors of 1 fac 0-8 and velocity increments of 5 ft/sec were used. For 

studying behavior along the positive real axis, particularly near the origin, starting locations 
along the real axis were used with initial velocities selected to capture any potential real, 
divergence root. Smaller relaxation factors and velocity increments/decrements were also 

required there. For the airfoil examples, the convergence criterion was |(s‘ -s„)| < 0.001 
rad/sec. Since A(s) is complex, the sets of eigenvalues, X(s l ), of A(s n ) for any converged 
root, s l = o' +ico l , are not in general a collection of real roots and complex conjugate root 
pairs. Checks confirmed that, for complex valued roots, their complex conjugates, o' -io) 1 , 
were also converged roots. Only the roots for the upper half s -plane are presented. 

The use of another shareware computational routine for calculation of the determinant of a 
general, complex matrix was very useful in gaining understanding of the location and 
development of singularities. Thus, by evaluating detjsl- ^(s)} along a ray in the s-plane, 
the singularities (e.g., poles and zeroes) of the system characteristic equation can be studied. 
Also, by evaluating det-f^I-^)} along a circle contour about a test point, the number of 

singularities within the circle may be studied by noting the number of 2 n radian phase 
changes. 

3.2 Wing model 


6 



The BAH wing model shown in Fig. 2 consists of the truncated normal mode equations 


+ 2 + mfii% = 


N 


p u2 1&,U i=l ’ N 

j = i 


( 12 ) 


where Q Uj are elements of the NxN complex GAF matrix, Q . The mode shapes, masses, and 
frequencies of all ten modes, A = 10, were extracted from the MSC/NASTRAN code. The 
viscous damping is zero for this case. Following transformation and identifying 
X T =(<£,,.. Eqs. 12 can be written in the same form as Eqs. 9 and 10 where 

B(s) = Oand K(s) = K s -\pU 2 Q(J,M) ■ The GAF matrix, Q, was obtained from a modified 
version of the DLAT module of the ISAC program (Interaction of Structures, Aerodynamics, 
and Controls) described by Adams and Hoadley 2 ’ 1 . The original DLAT module implements 
the DLM code described in Ref. 24, calculating the GAF matrix, Q(ik w ,M), for purely 
harmonic reduced frequencies. The modified DLAT module calculates the GAF matrix for 
generalized reduced frequency, s, in the manner described by Edwards 14 . Changes to the 
code are straightforward, with all occurrences of the real input reduced frequency variable, 
k in , replaced by the generalized complex variable, k = k in -iajj, signifying arbitrary I -plane 
motion with complex reduced frequency I = crjj + ik in . For example, in the DLM the kernel 
function of the integral equation to be solved involves functions such as Eq. A.26 in Ref. 24: 

/„(«„*,) = e ‘ k ' U ' £°[l-w(l + w 2 )" 1/2 y^clu (13) 


which is approximated by exponential series (Eq. A. 30 in Ref. 24) 


h Oh 9 ~ 


a„e 


-tnc +k; 


~{nc-ik x ) 


n = 1 


(14) 


where the real constants a„ and c are given in Ref. 24, k x = 6)r / U , and the variables U\ and r 
are lengths related to the surface integration of the kernel. Generalization is accomplished by 
generalizing the real variable, k v to the complex variable, 
= (a>-i<j)r /U = (k in -iojj)r/b. It is shown in Refs. 12, 13, and 15 that the stability of 

linear aeroelastic systems is determined by these generalized airloads. Portions of the system 
response which are not included (they are dependent upon initial conditions) do not affect 
system stability. 


In the Appendix of Ref. 16, Cunningham and Desmarais discuss the singularities of the 
subsonic kernel function, showing that Eq. 1 3 has no singularities in the finite part of the I - 
plane except for a logarithmic branch point at the origin. Similarly to the Theodorsen 
function, a branch cut along the negative real axis is required to make the function single- 
valued. Now, the approximation of Eq. 14 has zeroes at /c, =-inc and poles at /c, =±inc. 
Accounting for pole-zero cancellations, the approximation has poles located at 
s=-(blr}nc, n = 1,1 1 on the negative real axis in the s -plane. Ref. 16 discusses the 
effect these poles on the performance of several approximations like that of Eq. 14. As 
|arg(,s')- 7 T /2| increases towards nil, the approximations become unacceptable for 
determining accurate airloads. 
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4 RESULTS 


4.1 Airfoil models 

The root loci for the 2 DOF restrained airfoil model are shown in Fig. 3 for the c. g. at 37% 
chord and in Fig. 4 for the c. g. at 45% chord. The corresponding plots for the 3 DOF 
unrestrained airfoil model are shown in Figures 5 and 6. Table I summarizes the flutter 
speeds, t/^, and frequencies, co f , and the divergence speeds, U D (and frequencies, (0 D , for 
the 3 DOF cases) for these cases. For the 2 DOF case, divergence is aperiodic with the 
divergence root emerging from the origin at U D . For the 3 DOF case, the additional fuselage 
mass results in a low frequency oscillatory mode which becomes unstable at speeds near those 
of the restrained 2 DOF model divergence speeds. Rodden and Bellinger 5 refer to this mode 
as a ‘dynamic divergence’ mode, which seems appropriate. All of these root loci were 
calculated using the eigenvalue iteration method of Eq. 1 1 and the exact Theodorsen airloads 
for arbitrary complex values of J and thus are the exact roots of the flutter equation (within 
the error band due to the convergence criterion). All of the published calculations of these 
characteristics (see Refs. 2-6, 8, and 9) are close to those given in the table and differences are 
due to the various approximations to the airloads (c.f. Eqs. 2-5) used by these authors. Of 
course this is as it should be, as these stability characteristics are defined by the airloads for 
purely harmonic oscillations which is the foundation of all of these approximate methods. It 
should be expected that deviations between the exact root loci of Figs. 3-6 and those 
calculated by the approximate methods would become progressively larger as root locations, 
s, become further removed from the vicinity of the ia> axis. The ‘x’ symbols in Figs. 3-6 
give the coupled, wind-off root locations. (The virtual mass airloads are included.) The root 
locus branch in Figs. 3-6 originating from the coupled torsion mode near 25 rad/sec behaves 
similarly in each case. The mode initially becomes damped for increasing speed along with 
decreasing frequency due to interaction with the lower frequency plunge mode. The mode 
eventually becomes the flutter mode and achieves a maximum level of negative damping 
before becoming asymptotic to a complex frequency value near the ico axis. The root locus 
branch originating from the coupled plunge mode becomes monotonically more heavily 
damped in all cases. It never coalesces with its complex conjugate root to form real roots on 
the negative real axis. For both the 2 DOF and 3 DOF cases with U = 1000 ft I sec , the root 
is located at s = -100.87 + 130.89 rad/sec with the c.g. at 37% chord and at 
s = -1 13.65 + 136.97 rad/sec with the c.g. at 45% chord. 


The emergence of the divergence root as an additional, fifth root in the 2 DOF case was 
discovered by Edwards 12 . A simple explanation of this additional root is given by a low 
frequency model of a single degree of freedom pitching airfoil, which leads to the 
characteristic equation (c.f. Eq. 8-7 of Ref. 1) 


l-C(J) 


'uf 1 

u n 


c(s) = 0 with 


'stii 


(15) 


“A pole of this aeroelastic system occurs at values of s for which the coefficient in Eq. 16 is 
zero. Since C(s ) is purely real only on the positive real axis, poles can only occur there. 


Also, along the positive real axis, C(s ) decreases monotonically from a value of 1 .0 at |s| = 0. 
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Figure 3, Exact root loci as a function of 
airspeed, showing flutter and divergence: two 
degree of freedom airfoil, c. g. at 37% chord. 


Figure 5. Exact root loci as a function of 
airspeed, showing flutter and divergence: three 
degree of freedom airfoil, c. g. at 37% chord. 
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Figure 4. Exact root loci as a function of 
airspeed, showing flutter and divergence: two 
degree of freedom airfoil, c. g. at 45% chord. 



Figure 6. Exact root loci as a function of 
airspeed, showing flutter and divergence: three 
degree of freedom airfoil, c. g. at 45% chord. 


to 0.5 at \s\ = °° . Hence, a pole cannot occur for U <U D and for U >U D , only one real pole 

can occur. This mode produces the motion of the diverging airfoil and occurs in addition to 
the 2n structural poles.” 12 Another explanation of the possibility of additional roots is based 
on a fundamental theorem of algebra: an nth order constant coefficient polynomial equation 
(characteristic equation) has n and only n roots, but if the coefficients are not constant, as in 
Eqs. 1 and 9, no such guarantee exists. 


The emergence of the divergence root at the origin for the 2 DOF case with the c.g. at 45% 
chord is illustrated in Fig. 7. The complex determinant of the stability matrix is evaluated 
along rays in the 5 -plane just above and just below the real axis at speeds bracketing the 
divergence speed, U D = 216.6 ft/sec. The development and emergence of the divergence pole 
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is clearly seen. The effect of the branch 
cut along the negative real axis is also 
evident in the jump in phase there. Figure 
8 shows the complex stability determinant 
at a higher speed, U = 315 ft/sec, with the 
divergence pole located at s = +1 1 rad/sec. 


The ability of aeroelastic models utilizing 
state-augmentation to accurately predict 
static divergence speeds has become 
accepted. Models akin to that of Eq. 5, 
termed Rational Function 

Approximations 25 (RFA), generally 
employ a sequence of negative, real 
eigenvalues (roots) bracketing the 
frequency range of interest. The resulting 
additional roots have been termed ‘aerody- 
namic lag roots’ 7 . They tend to cluster 
about the negative real axis and as speed increases towards u D the root closest to the origin 
migrates to the right, crossing onto the positive real axis at U D . Edwards 12 illustrates this for 
an airfoil model with stability behavior similar to that of the above 2 DOF model (with the 
c.g. at 37% chord). 

The 3 DOF airfoil case of Figs. 5 and 6 presents a different situation regarding the number of 
singularities in the problem. Here the stability matrix (Eqs. 6 and 9) is 6x6 and at least six 
roots would be expected. However, the complex determinant shown in Fig. 9 indicates the 
presence of a seventh root at the low speed of U = 5 ft/sec for this case with the c.g. at 45% 
chord. An isolated singularity (pole) is evident at the origin, along with a pair of damped, 
complex conjugate roots which constitute the dynamic divergence mode at U = 232.9 ft/sec 
(see Fig. 6). Fig. 10 shows the presence of three real roots at the higher speed of 
(7 = 315 ft/sec; the additional pole at the origin and the two poles resulting from the merger of 
the dynamic divergence mode poles on the positive real axis at U = 308. 15 ft/sec (see Fig. 6). 
(The four other roots of the bending and torsion modes are not displayed in Figs. 9 and 10.) 
Comparison of Figs. 8 and 10, at the same speed, is instructive in illustrating the differences 
between the 2 DOF and 3 DOF airfoil cases. Physically, the seventh root at the origin 
accounts for a quasi-static airfoil motion with constant plunge velocity for this unrestrained 
airfoil case. 

Further calculations, not shown, explored the transition from the 3 DOF model to the 2 DOF 
model by calculating root loci for increasing values of the fuselage mass, m f , in Eq. 6. As 

m f increases, the bending and torsion mode loci of the 3 DOF models of Figs. 5 and 6 

approach those of the 2 DOF models of Fig. 3 and 4. Also, the ‘oval’ locus of the dynamic 
divergence mode near the origin becomes smaller, shrinking to the origin in the limit of 
infinite m f , and the dynamic divergence speed approaches the static divergence speed of the 

2 DOF model. Rodden 26 also notes this behavior. Of course, in the limit of infinite m f , the 

3 DOF model with its 7 singularities becomes the 2 DOF model with its 5 singularities. 



2 DOF 

3 DOF 

C.g., % c 

37 

45 

37 

45 

U f , ft/sec 

257.1 

169.1 

284.1 

159.5 

(O f , rad/sec 

15.64 

16.07 

16.84 

17.37 

U D , ft/sec 

217.0 

216.6 

232.9 

215.2 

0 ) D , rad/sec 

- 

- 

7.29 

7.30 


Table I. Flutter and divergence 
characteristics of restrained and 
unrestrained airfoil models. 
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Figure 7. Determinant of stability matrix versus 
Real(s): two degree of freedom airfoil, c.g. at 
45% chord. 




Figure 9. Determinant of stability matrix versus 
Real(s): three degree of freedom airfoil and 
fuselage, c.g. at 45% chord, U = 5 ft/sec. 



Figure 8. Determinant of stability matrix versus 
Real(s): two degree of freedom airfoil, c.g. at 
45% chord, U = 315 ft/sec. 




Figure 10. Determinant of stability matrix 
versus Real(s): three degree of freedom airfoil 
and fuselage, c.g. at 45% chord, U = 315 ft/sec. 


4.2 BAH wing model 

For three-dimensional wing configurations, like the BAH wing, the stability matrix in Eq. 9 
again contains non-constant elements, i. e. the airloads, Q(s,M). The same procedure 
described above for the airfoil cases was used to determine the root loci of the model. Ten 
modes were included resulting in a 20x20 general, complex stability matrix. The wind-off, 
coupled structural mode locations were used as starting values for the oscillatory roots. The 
velocity, for this M = 0 case, was started at 100 ft/sec and incremented by 2 percent following 
convergence at each speed. The reference sea-level density is p ref - 1.1468xl0' 7 sl(in)/in. 3 

In this case, the virtual mass in not included in the wind-off calculations, accounting for the 
small jump from the ‘x’s at the start of each locus. The iterative relaxation method of Eq. 1 1 
was again very robust for the oscillatory modes and a relaxation factor of 0.80 was again 


11 



used. For the real divergence mode root, speed increments or decrements of 0.5 percent and a 
relaxation factor of 0.02 were required. The convergence criterion was 

|(s' -s„)| < 0.03 rad/sec. 

The root loci for the first seven modes and the divergence mode are shown in Fig. 11. The 
lowest flutter speed involves the second mode with = 1059.7 ft/sec and 6^ =19.39 

rad/sec. A real, non-oscillatory root was found for 1 450 < U <2036 ft/sec with divergence 
occurring at U D = 1647.25 ft/sec. Flutter modes are also found at higher speeds for modes 3, 
5, and 6 at U =2449, 2217(5280), and 6000 ft/sec respectively. For mode 5, looping of the 
locus produces two flutter speeds in the speed range covered. Such looping of the stability 
root loci is reminiscent of looping of airfoil load coefficients for lightly damped oscillations at 
high reduced frequencies shown by Cunningham and Desmarais 16 . Fig. 12 provides another 
viewpoint of the flutter instability from the determinant of the stability matrix versus ia> at 
U = 1056 ft/sec, near the flutter speed. The presence of the flutter mode is clearly evident at 
io) = ia f . 

This example with its very coarse structural and aerodynamic paneling cannot be expected to 
yield high accuracy aeroelastic results, especially for the higher frequency modes and high 
values of reduced frequency. It is, however, a very good tutorial checkcase and has been well 
used to highlight differences obtained from alternative analyses. Fig. 1 1 illustrates that the 
‘exact’ aeroelastic root loci for the oscillatory modes of this problem are well behaved, with 
no discontinuities, jumps, or bifurcations. (The present results are termed ‘exact’ in a sense 
that is discussed below.) Crossings of the loci do occur, but for very different speeds and lead 
to no ambiguity in interpreting the results, in contrast to the difficulties encountered in the k 
and p-k solution methods. 21 It is the opinion of the authors that the s - plane presentation, 
with velocity as a parameter, gives a much clearer view of system stability than U-g and 
U-f diagrams. There is also no need to introduce decay rates as damping values for non- 
oscillatory real roots in order to interpret the plotted results. 

An important result of this study is the insight gained of the behavior of the root loci along the 
real axis, particularly for the negative real axis. The loci of the first three oscillatory modes in 
Fig. 1 1 are shown approaching the negative real axis. Numerical convergence difficulties 
prevented further calculations for these loci and also for the locus of the real root for speeds 
U < 1450 ft/sec. To understand the numerical issues involved, the determinant of the stability 
matrix versus Re(s) = o was calculated for several velocities. Figs. 13 and 14 show the 

determinant for [7 = 1450 and U = 1647.25 (U D ) ft/sec, respectively. The source of the 

numerical difficulty is immediately apparent in the numerous singularities (poles and zeroes) 
occurring along the negative real axis starting about -7 rad/sec. For U = 1450 ft/sec. Fig. 13, 
an isolated pole has migrated out of the numerically difficult region, becoming the divergence 
mode, at U D = 1647.25 ft/sec. Fig. 14. 

We return now to discuss the intent in referring to the current stability loci results as ‘exact’ 
results. The GAF matrix, Q(s,M ), produced by the generalized DLAT module (the DLM 
code) contains the exact solutions for the linear matrix equation resulting from the several 
approximations utilized in the DLM formulation. The subsequent use of Q(s,M) in 
determining the root loci of the stability matrix of Eq. 9 throughout the s -plane is exact (a 
small error is due to the convergence criterion of the eigenvalue iteration method). Thus, any 
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Figure 11. ‘Exact’ root loci as a function of 
airspeed showing flutter and divergence: 
transport (BAH) wing model, M = 0, 



lmag(s), rad/sec 


Figure 12 . Determinant of stability matrix 
versus Imag(s): transport (BAH) wing model, U 
= 1056 ft/sec. 
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Figure 13. Determinant of stability matrix 

versus Real(s): transport (BAH) wing model, U Figure 14. Determinant of stability matrix 

= 1450 ft/sec. versus Real(s): transport (BAH) wing model, U 

= 1647.25 ft/sec. 


error involved in the current root loci calculations is that introduced by the reduction of the 
linear potential equation boundary value problem (BVP) to a linear matrix equation. The 
solution of the BVP is expressed as a singular integral equation 11,16 involving a kernel 
function whose singular behavior has been discussed above 16 . Two main methods of solution 
of this integral equation are the Kernel Function Method (see, e.g. Ref. 16) and the Doublet 
Lattice Method (Ref. 1 1). Both methods reduce the singular integral equation problem to a 
linear algebraic problem for the pressure loads on the surface and both methods use several 
approximations to accomplish this reduction. Both utilize exponential series similar to Eq. 14 
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to approximate expressions like Eq. 13. The effect of this upon the singular behavior of the 
problem is to replace the logarithmic branch point singularity of the kernel function with 
isolated singularities (poles, or lag roots) typically located along the negative real axis as 
discussed above. Other approximations employed in the DLM are 
the discretization into the lattice of doublet lines and a parabolic approximation of the kernel 
along doublet lines. Approximations used in the Kernel Function Method are the truncation of 
the series approximation of the pressure distribution and the choice of the pressure loading 
functions to be used in the truncated series. It is with this understanding of the 
approximations introduced by the generalized DLM solution that the root loci computed by 
the current eigenvalue iteration method are termed ‘exact’. 

The similarity between RFA models and the exponential series approximation, Eq. 14, used in 
the DLM suggests that the singularities in Figs. 13 and 14 for Re(.S’) < -6 rad/sec derive from 
this series approximation. (The RFA given by the transform of Eq. 5 has negative, real poles 
as does Eq. 14.) Clusters of lag roots on and near the negative real axis similar to those seen 
in Figs. 13 and 14 are a well known feature in RFA flutter analyses. 

These discussions of DLM solutions in the vicinity of the negative real axis lead naturally to 
the issue 4 of the origin of the divergence mode in Fig. 11. Does the divergence mode derive 
from the bending mode or from a lag root? From Fig. 11 it is apparent that the bending mode 
and its complex conjugate root do merge onto the real axis at a speed slightly higher than 
1392 ft/sec and that a real root does emerge from the cluster of roots at a speed slightly below 
1450 ft/sec. The change in the nature of the singular behavior of the kernel function resulting 
from the use of exponential approximations like Eq. 14 must also be acknowledged here. 
That is, the complication of the dense collection of poles and zeroes along the negative real 
axis is an artifact of the solution approximations, masking the behavior of the true exact 
solution in this region. Along with the above discussion of the inaccuracy of solutions for 
these highly damped values of complex frequency, this would appear to render the issue of 
the origin of the divergence root moot. In any case, the root has not been tracked for 
1392 <U< 1450 ft/sec. 

The robustness of the GAAM eigenvalue iteration procedure in tracking the loci of the 
oscillatory structural modes of Fig. 11 was investigated. At issue is the maximum increment 
in velocity allowable in order to track the proper locus. For the two lowest frequency modes 
near 12 and 22 rad/sec all velocity increments up to 1000 ft/sec were able to track the loci 
over the range of velocity shown in the figure. For the third, fourth, and fifth modes near 46, 
74, and 94 rad/sec the maxima were 500, 200, and 100 ft/sec respectively. Of course these 
maximum increments are related to the closeness of the roots on the loci for corresponding 
velocities; c. f. the fourth and fifth mode loci near s ~-10 + ?60 rad/sec for U -4500 ft/sec. 

Finally, the current Generalized Aeroelastic Analysis Method using the generalized DLM 
code is not restricted to the incompressible, M = 0, condition of the BAH wing example. 
The main application of the DLM is in fact to lifting surface configurations in subsonic flow. 
To this end, root loci of the BAH wing model flutter mode versus Mach number were 
calculated. For this example, the sound speed is assumed to bea„ = 1500 ft/sec and U = Ma„ . 
Fig. 15 shows the Flutter mode locus for two densities: the solid line is for the density of the 
above example and the dashed line is for a 50 percent higher density. The flutter Mach 
number (and velocity) decreases from 0.778 to 0.662 with increasing density. The symbols 
denote 0.10 Mach number increments. Note that each point on these loci are ‘exact’ matched 
point values of frequency and damping and are directly calculated by GAAM with no root- 
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sorting, complex mode tracking logic, or 
dataset interpolation required. 

5 DISCUSSION 

All of these ( M = 0) cases have been 
studied extensively by a number of ® 

researchers using various assumptions in | 

the use of harmonic oscillation airloads for 2 

the several analysis methods described | 

above. The 2 DOF airfoil cases have been 
analyzed in Refs. 2-4 and the 3 DOF cases 
have been analyzed in Refs. 3, 5, 6, and 9. 

The BAH wing model has been studied in 
Refs. 3, 7, and 8. There are many 
differences in solution details between the 
exact results given in Figs. 3-6, the ‘exact’ Figure 15. ‘Exact’ root loci of flutter mode as a 

results given in Fig. 11, and those reported function of Mach number and density: transport 

for these approximate solutions. Rather (BAH) wing model. 

than focusing on each detail, the discussion to follow offers comments on the origin of these 
differences. All solutions produce good agreement for the flutter and divergence properties of 
the cases. This to be expected as all methods are based upon airloads which are accurate for 
purely harmonic oscillations (flutter and dynamic divergence) and for steady flow (static 
divergence). Differences in these properties are due solely to the differences utilized in these 
references for interpolating or approximating the oscillatory airloads. For speeds which are 
well removed from these flutter and divergence speeds and for highly damped or undamped 
modes, there are curious behaviors reported for the aeroelastic root loci (as shown in 
U-g and U-f plots). These include ‘discontinuous roots’ 2 ’ 3 , speeds for the ‘activation of 
aerodynamic lag roots’ 3,8 , ‘bifurcations’ of the roots’’ 4 , inconsistency over the origin of 
certain of the roots 4 , and oscillatory roots ceasing to exist above critical speeds 4 . 

These behaviors can be attributed to the use of harmonic oscillation airloads to infer loads for 
complex frequencies well removed from the ico axis. Edwards 12 ' 15 proved that airloads 
derived for harmonic oscillation conditions can be used to derive airloads off the i(0 axis by 
appealing to analytic continuation. There are, however, practical limits to this process. 
Regarding the airfoil cases, note that all of the exact loci shown in Figs. 3-6 are continuous 
and have no discontinuities or jumps. No lag roots become ‘active’ at any speed and the only 
type of root ‘bifurcation’ which occurs is for the dynamic divergence mode roots joining the 
real axis in Figs. 5 and 6. All of these types of features shown in Refs. 2-6, and 16 which 
differ from the corresponding root loci in Figs. 3-6 are spurious and are introduced by the 
airload modeling approximations used or by the aeroelastic solution methods. 

A clarification of terminology is needed regarding use of the term ‘bifurcation’ to describe the 
merging of the dynamic divergence mode roots in Figs. 5 and 6 onto the positive real axis 
near s~ 6 -8 rad/sec. In nonlinear system theory, locations in parameter space at which a 
qualitative change in response occurs are termed bifurcation points. For the linear aeroelastic 
systems considered in this paper no such qualitative change in response occurs. These 
locations should be termed break-in (or breakaway) points. 
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The issue of an additional real root emerging from the origin for restrained airfoil cases, Eq. 
16, has also been addressed by Heeg 27 who analyzed a 1 DOF pitching airfoil model using an 
unsteady discrete vortex lattice method and eigenmode analysis. By varying the mass ratio, 
elastic axis, and radius of gyration, cases are shown in which the divergence mode is an 
additional eigenmode (as in Figs. 3 and 4) and in which the divergence mode results from the 
pitching mode root locus joining the real negative axis at a break-in point with subsequent 
divergent behavior. In the latter case, there would be no additional singularity. Ref. 27 also 
provides experimental verification of the former case from low speed wind tunnel model tests. 

Regarding the BAH wing case (and by extension, three-dimensional lifting surface 
applications in general), the number of oscillatory frequencies for which the loads are 
available is limited in practice and can be insufficient to accurately define loads well removed 
from the id) axis. Also, Ref. 16 shows that sin-gularities introduced by the exponential series, 
such as Eq. 14, in the kernel function solution lead to errors rendering the resulting ‘exact’ 
solutions unusable for significant regions of the s -plane. For a 12-term series 16 solutions for 
|arg(T)|</r/4 or |arg(s)| > 3/r/4 become inaccurate or uneconomical, or both. An 11 -term 

series 24 commonly used in DLM codes (and in the code used in the present study) was 
evaluated and found to be acceptable for arg(T) = 5^/8 but unacceptable for arg(T) = Zn / 4 . 
Thus, for the BAH wing all root loci in the vicinity of the negative real axis which are 
presented in Refs. 3, and 7-9 or shown in Fig. 11 have questionable accuracy. 

Most of the curious behaviors of aeroelastic root loci mentioned above occur for heavily 
damped/undamped or real values of .S' where this issue of accuracy of the airloads is most 
pronounced. On the other hand, the root loci of Fig. 11 obtained with the generalized DLM 
are the exact solutions to the linearized algebraic problem resulting from the approximations 
to the kernel functions discussed above. The harmonic oscillation DLM airloads used in Refs. 
3, and 7-9 are all based upon this same linearized algebraic problem. Thus the ‘exact’ root 
loci of Fig. 1 1 can be used to evaluate corresponding solutions from these references. Fig. 1 1 
shows that the ‘exact’ root loci are continuous and have no discontinuities or jumps. No 
bifurcations of roots occur and only one real root participates in the static divergence of the 
wing. Any differences among the ‘exact’ root loci of Fig. 1 1 and those from Refs. 3, and 7-9 
are introduced by the airload modeling approximations or by the aeroelastic solution methods 
used in these references. 

Desmarais and Rowe 2s present alternative kernel approximation algorithms which are tailored 
to enable accurate approximations throughout the several regions of the complex plane. In 
place of exponential series like Eq. 14, the algorithms utilize Neuman series, continued 
fractions, asymptotic series, and Modified Functions. To the authors knowledge, they have 
not been implemented in an aeroelastic analysis program. 

It is interesting to contrast this predicted behavior of the bending mode of the BAH wing 
model approaching divergence with that of the forward swept wing wind tunnel model 
described in Ref. 29. There, one of three cantilevered elastic models was lost due to static 
divergence while testing at M = 1.05 and at 88 percent of the predicted divergence dynamic 
pressure. Significantly, the bending strain gage signal clearly indicated the presence of the 2 
Hz. bending mode throughout the event. It is possible that this event is similar to that shown 
in Fig. 3 where divergence occurs due to the emergence of an additional singularity at the 
origin while the bending mode remains distinct and oscillatory. 
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Two main points resulting from the present investigation are: (i) the possibility of additional 
singularities occurring in aeroelastic analyses over and above those introduced by the 
structural dynamics and (ii) inherent limitations in computational codes commonly used for 
aeroelastic analyses for the study of behavior well removed from the ico axis. It must be 
emphasized that only a small number of examples have been studied. No attempt has been 
made to delineate boundaries of parameters which exhibit the novel singular behavior shown 
in Figs. 3 and 4. The studies of Ref. 27 imply that sharp boundaries may not be found. 
Problems with the accuracy of the airloads for heavily damped wing motions notwithstanding, 
the BAH wing example appears to show that divergence for this case does occur subsequently 
to the bending mode frequency dropping to zero. It has been past practice to assume that this 
pattern must be involved in cantilevered wing divergence 29 . Results of this study (and the 
experimental results of Refs. 27 and 29) show that this assumption must be called into 
question. There is also the possibility of new insights which may be beneficial for subcritical 
divergence testing methods. 

The following commentary provides further insight on issues dealt with herein: “An incorrect 
solution to an ‘unrestrained divergence problem’ was given in Ref. 10 using the static method 
and has been corrected in an errata 30 . In the errata, it was observed that in Ref. 10 the 
presumed instability occurred at a dynamic pressure at which the mean structural axis of the 
vehicle remained aligned with the freestream velocity vector rather than at the dynamic 
pressure of an actual physical structural instability. The incorrect quasi-static solution was 
applied to two example airfoils with an additional fuselage degree of freedom 5 (c.f. the cases 
treated in Figures 5 and 6 herein) and compared with a dynamic solution. In one case the 
quasi-static divergence speed was lower by 1.7% than the dynamic solution; in the second 
case the difference was lower by 6.8%. Since the errata had not yet been published, the 
explanation for the discrepancies offered 5 was also in-correct. The conclusion of the errata 
was that ‘divergence of an unrestrained vehicle should always be investigated by dynamic 
stability methods.’” (private communication, W. P. Rodden, July 2002) Dykman and 
RoddeiT 1 also treat the transient response of an unrestrained, flexible vehicle and compare the 
dynamic ‘correct’ solution with those obtained using modal residualization, modal truncation, 
and quasi-steady aerodynamics. 

It must be emphasized that the use of the Generalized Aeroelastic Analysis Method for the 
calculation of aeroelastic root loci as illustrated in the above examples has not cast doubt on 
other methods of determining flutter and divergence characteristics (velocity and frequency). 
However, the ability of performing aero-elastic analysis without the necessity of root-sorting 9 , 
root-searching 4 , lining-up 9 , or reduced-frequency-sweep’’ techniques is attractive. It is 
important to understand that each converged root locus value shown in Figs. 3-6, 11, and 15 
is a ‘matched-point’ solution. By appropriate specification of density, Mach number, and 
speed direct calculations of such ‘matched-point’ aeroelastic root loci can be computed for 
constant altitude with varying Mach number and velocity, or for constant Mach number with 
varying altitude and velocity. 

6 CONCLUSIONS 

The Generalized Aeroelastic Analysis Method (GAAM) is applied to the analysis of three 
well-studied checkcases: restrained and unrestrained airfoil models, and a wing model. An 
eigenvalue iteration procedure is used for converging upon roots of the complex stability 
matrix. For the airfoil models, exact root loci are given which clearly illustrate the nature of 
the flutter and divergence instabilities. The singularities involved are enumerated, including 
an additional pole at the origin for the unrestrained airfoil case and the emergence of an 
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additional pole on the positive real axis at the divergence speed for the restrained airfoil case. 
Inconsistencies and differences among published aeroelastic root loci and the new, exact 
results are discussed and resolved. The generalization of a Doublet Lattice Method computer 
code is described and the code is applied to the calculation of root loci for the wing model for 
incompressible and for subsonic flow conditions. The error introduced in the reduction of the 
singular integral equation underlying the unsteady lifting surface theory to a linear algebraic 
equation is discussed. Acknowledging this inherent error, the solutions of the algebraic 
equation by GAAM are termed ‘exact.’ The singularities of the problem are discussed and 
exponential series approximations used in the evaluation of the kernel function shown to 
introduce a dense collection of poles and zeroes on the negative real axis. Again, 
inconsistencies and differences among published aeroelastic root loci and the new, ‘exact’ 
results are discussed and resolved. In all cases, aeroelastic flutter and divergence speeds and 
frequencies are in good agreement with published results. The GAAM solution procedure 
allows complete control over Mach number, velocity, density, and complex frequency, thus 
all points on the computed root loci can be matched-point, consistent solutions without 
recourse to complex mode tracking logic or dataset interpolation, as in the k and p-k solution 
methods. 

7 ACKNOWLEDGEMENT 

The authors wish to acknowledge Dr. William P. Rodden of La Canada Flintridge, California 
for his advice and discussions on the topics treated in the paper over the past quarter-century 
and for the commentary included above. The Aeroelasticity community has been well served 
in his development and/or advancement of these computational checkcases. 

8 REFERENCES 

[1] Bisplinghoff, R. L., Ashley, H., and Halfman, R. L., Aeroelasticity, Addison-Wesley, 
Reading, Mass., 1955. 

[2] Rodden, W. P. and Bellinger, E. D., “Aerodynamic Lag Functions, Divergence, and the 
British Flutter Method,” Journal of Aircraft, Vol. 19, No. 7, 1982, pp. 596-598. 

[3] Chen, P. C., “A Damping Perturbation Method for Flutter Solution: The g-Method,” 
Proceedings of the International Forum on Aeroelasticity and Structural Dynamics, Pt. 1, 
CP 1 999-209 136/PT1, NASA, 1999, pp. 403-413. 

[4] van Zyl, L. H., “Aeroelastic Divergence and Aerodynamic Lag Roots,” Journal of 
Aircraft, Vol . 38, No. 3, 2001, pp. 586-588. 

[5] Rodden, W. P. and Bellinger, E. D., “Unrestrained Aeroelastic Divergence in a Dynamic 
Stability Analysis,” Journal of Aircraft, Vol. 19, No. 9, 1982, pp. 796-797. 

[6] van Zyl, L. H., “Unrestrained Aeroelastic Divergence and the p-k Flutter Equation,” 
Journal of Aircraft, Vol . 38, No. 3, 2001, pp. 588-590. 

[7] Rodden , W. P., and Stahl, B. “A Strip Method for Prediction of Damping in Subsonic 
Wind Tunnel and Flight Flutter Tests,” Journal of Aircraft, Vol. 6, Jan.-Feb. 1969, pp. 9-17. 

[8] van Zyl, L. H., “Divergence and the Flutter Equation,” Proceedings of the International 
Forum on Aeroelasticity and Structural Dynamics, Asociacion de Ingenieros Aeronauticos de 
Espana, Vol II, 2001, pp. 401-412. 

[9] Rodden, W. P., and Johnson, E. H., ‘MSC/NASTRAN Version 68 Aeroelastic Analysis 
User’s Guide, MacNeal-Schwendler Corp., Los Angeles, CA, 1994. 

[10] Rodden, W. P., “Aeroelastic Divergence of Unrestrained Vehicles,” Journal of Aircraft, 
Vol. 18, No. 12, 1981, pp. 1072-1073. 


18 



[11] Albano, E., and Rodden, W. P., “A Doublet-Lattice Method for Calculating Lift 
Distributions on Oscillating Surfaces in Subsonic Flows,” AIAA Journal, Vol. 7, No. 2, 1969, 
pp. 279-285. 

[12] Edwards, J. W., “Unsteady Aerodynamic Modeling and Active Aeroelastic Control,” 
Stanford University SUDAAR 504, Feb. 1977. 

[13] Edwards, J. W., “Unsteady Aerodynamic Modeling for Arbitrary Motions,” AIAA 
Journal, Vol. 15, Vo. 4, April 1977, pp. 593-595. 

[14] Edwards, J. W., “Application of Laplace Transform Methods to Airfoil Motion and 
Stability Calculations,” AIAA Paper No. 79-0772, St. Louis, MO, April 1979. 

[15] Edwards, J. W., Ashley, H., and Breakwell, J. V., “Unsteady Aerodynamic Modeling for 
Arbitrary Motions,” AIAA Journal, Vol. 17, No. 4, April 1979, pp. 365-374. 

[16] Cunningham, H. J. and Desmarais, R. N., “Generalization of the Subsonic Kernel 
Function in the s-Plane, With Applications to Flutter Analysis,” NASA TP-2292, March 
1984. 

[17] Cunningham, H. J., “Analysis of Preflutter and Postflutter Characteristics With Motion- 
Matched Aerodynamic Forces,” NASA TP 1232, July 1978. 

[18] Hassig, H. J., “An Approximate True Damping Solution of the Flutter Equation by 
Determinant Iteration,” Journal of Aircraft, Vol. 8, No. 1 1, 1971, pp. 885-889. 

[19] Rodden, W. P., Harder, R. L., and Bellinger, E. D., “Aeroelastic Addition to 
NASTRAN,” NASA CR 3094, 1 979. 

[20] Jones, W. P., “Aerodynamic Forces on Wings in Non-Uniform Motion,” R and M 2117, 
1945, Aeronautical Research Council. 

[21] Eldred, M. S., Venkayya, V. B., and Anderson, W. J., “Mode Tracking Issues in 
Aeroelastic Analysis,” AIAA Journal, Vol. 33, No. 7, 1995. 

[22] Goodman, C., “Accurate Subcritical Damping Solution of Flutter Equation Using 
Piecewise Aerodynamic Function,” Journal of Aircraft, Vol. 38, No. 4, July-August 2001, pp. 
755-763. 

[23] Adams, W. M. and Hoadley, S. T., “ISAC, A Tool For Aeroservoelastic Modeling and 
Analysis,” NASA TM 1 0903 1 , December 1993. 

[24] Giesing, J. P., Kalman, T. P., and Rodden, W. P., “Subsonic Unsteady Aerodynamics for 
General Configurations. Part I, Volume I - Direct Application of the Nonplanar Doublet- 
Lattice Method,” AFFDL-TR-71-5, Pt. I, Vol, I, U. S. Air Force, November 1971. 

[25] Karpel, M. and Hoadley, S. T., “Physically Weighted Approximations of Unsteady 
Aerodynamic Forces Using the Minimum State Method,” NASA TP 3025, March 1991. 

[26] Rodden, W. P, Secondary Considerations of Static Aeroelastic Effects on High- 
Performance Aircraft, Paper No.3 in Static Aeroelasticity in Combat Aircraft, AGARD 
Report No. 725, 1986. 

[27] Heeg, J., “Dynamic Investigation of Static Divergence: Analysis and Testing,” Ph. D. 
Dissertation, Duke University, Durham, North Carolina, May 2000. 

[28] Desmarais, R. N. and Rowe, W. S, A Method for Computing the Kernel of the 
Downwash Integral Equation for Arbitrary Complex Frequencies, AIAA Paper No. 84-0983, 
May 1984. 

[29] Ellis, J. W. Dobbs, S. K, and Miller, G. D., “Structural Design and Wind Tunnel Testing 
of a Forward Swept Fighter Wing,” AFWAL-TR-80-3073, July 1980. 

[30] Rodden, W. P., “Aeroelastic Divergence of Unrestrained Vehicles,” Journal of Aircraft, 
Vol. 21, No. 1, January 1984, pp. 94-96. 

[31] Dykman, J. R. and Rodden, W. P., Structural Dynamics and Quasistatic Aeroelastic 
Equations of Motion, Journal of Aircraft, Vol 37, No. 3, May-June 2000, pp. 538-542. 


19 



